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We show that backflow correlations in the variational wave function for the Hubbard model greatly 
improve the previous results given by the Slater- Jastrow state, usually considered in this context. 
We provide evidence that, within this approach, it is possible to have a satisfactory connection with 
the strong-coupling regime. Moreover, we show that, for the Hubbard model on the lattice, backflow 
correlations are essentially short range, inducing an effective attraction between empty (holons) and 
doubly occupied sites (doublons). In presence of frustration, we report the evidence that the metal 
to Mott-insulator transition is marked by a discontinuity of the double occupancy, together with a 
similar discontinuity of the kinetic term that does not change the number of holons and doublons, 
while the other kinetic terms are continuous across the transition. Finally, we show the estimation 
of the charge gap, obtained by particle-hole excitations a la Feynman over the ground-state wave 
function. 

PACS numbers: 71.10.Fd, 71.30.+h 



I. INTRODUCTION 

Backflow correlations were introduced by Feynman 
and Coheni in order to obtain an accurate description 
of the excitation spectrum of liquid He''. Indeed, by ac- 
counting for these correlations, it was possible to have a 
marked improvement in the one-phonon dispersion curve 
and obtain an accurate description of the roton excita- 
tions. Indeed, it was realized that a picture of inde- 
pendent elementary excitations had to be extended to 
account for their interactions. Essentially, the current 
associated with a Feynman excitation involves a contri- 
bution from the backflow of the fluid around it. Feynman 
and Cohen showed that, in the simplest case of a particle 
tearing through the liquid at a given velocity, the pattern 
of induced longitudinal current far away from the parti- 
cle has a dipolar formii Then, the concept of backflow 
has been extended to weakly-correlated electron systems, 
where it turned out to be crucial for improving the de- 
scription of the electron gas in two and three dimensions, 
in particular, for correlation energies and pair distribu- 
tion functions Further applications concerned metallic 
Hydrogen^ and small atoms or molecules where signif- 
icant improvements in the total energy have been ob- 
tained. 

Very recently, backflow correlations have been success- 
fully applied to strongly-correlated electron systems, such 
as the Hubbard model and its generalizations to include 
frustrating terms, improving the variational wave func- 
tions used so far to approximate the exact ground state. 
In particular, in presence of frustration (e.g., for the t—t' 
Hubbard model on the square^ and the triangular lat- 
tices^) a spin-liquid phase has been stabilized only thanks 
to backflow correlations. This new term allowed us also 
to observe a renormalization of the underlying Fermi sur- 



face to perfect nesting at the metal-insulator transition.— 
Backflow correlations have been also recently applied to 
the infinite-t/ Hubbard model to study the stability of 
the Nagaoka ferromagnetism.*^ 

When compared to former variational wave functions 
for frustrated Hubbard modelsjS, (short-range) backflow 
correlations significantly improved the description of 
ground-state properties with respect to density-density 
Jastrow factors^ and holon-doublon binding factorsJ^ 
However, also in presence of backflow terms, the long- 
range Jastrow factor is still a fundamental ingredient of 
the wave function in order to describe the metal-insulator 
transition occurring at a finite value of the ratio between 
the electron-electron repulsion U and the hopping inte- 
gral t. 

In this paper, we show the accuracy of backflow corre- 
lations and compare them with the ^-matrix (i.e., strong- 
coupling) expansion] ^^1^^ The latter approach becomes 
exact for U/t ^ oo and define an iterative procedure 
which results in an expansion in powers of t/U. As a 
consequence, the (exact) wave function will be an ad- 
mixture of electronic configurations with weights pro- 
portional to (t/U)", where n is the number of doubly- 
occupied sites. Unfortunately, dealing with such a state 
is a difficult task, and, at present, only few attempts 
have been pursued in this directioufi^ On the contrary, 
our approach based upon backfiow correlations define a 
many-body wave function that can be easily treated by 
standard Monte Carlo methods. Most importantly, it 
remains accurate also for intermediate electron-electron 
repulsion, allowing a precise description of the physics 
down to the metal-insulator transition. ^'^ Here, we show 
that the energy gain due to backflow correlations in the 
insulating phase is mainly related to a better description 
of the local hopping processes changing the total num- 
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bers of holon-doublon couples. We discuss the exten- 
sion of backflow correlations to long-range distances and 
how the behavior of the Jastrow factor is affected by the 
presence of backflow correlations. We present the results 
both for one and two dimensions, with the variational 
wave function containing no magnetic terms. Finally, we 
discuss also how the variational method allows us to cal- 
culate the charge gap in the insulating phase and present 
results for the one-dimensional case. 

The paper is organized as follows: in Sec. HIl we in- 
troduce the Hamiltonian and we describe our variational 
wave function; in Sec. Illli we give a short review of the 
strong-coupling approach based upon the S matrix, to 
give a justification of the backflow terms; in Sec. IIVI we 
describe the backflow correlations; in Sec. |Vl we present 
some results for the metal-insulator transition and the 
insulating phase; finally, in Sec. I VII we draw our conclu- 
sions. 



II. MODEL AND VARIATIONAL WAVE 
FUNCTION 

We consider the Hubbard model with extended hop- 
ping in one-dimensional (ID) and two-dimensional (2D) 
square lattices. 

^ = - X! *y ^.^S,'^ + ^Y1 "^f^^,;: (1) 

i,j,(T i 

where ^ {^i denotes the creation (annihilation) op- 
erator of one electron on site i with spin a 4, rii^^ = 

A a'^i a is the electron density, tij the hopping parame- 
ters, and U the on-site Coulomb repulsion. In the follow- 
ing, we will consider nearest- and next-nearest-neighbor 
hoppings that will be respectively denoted by t and t' in 
2D and by ti and ^2 in ID. We will focus our attention 
on the half-filled case with L electrons on L sites. 

Both metallic and insulating phases can be constructed 
in a variational approach. In a first step, one constructs 
uncorrelated wave functions given by the ground state 
|BCS) of a superconducting Bardeen-Cooper-Schrieffer 
(BCS) Hamihonian^ii^ 

Hbcs = ^qcl^aCq,^ + ^g^.t^-g^ + ^-"^-^ 

q,a q 

where both the free-band dispersion tq and the pairing 
amplitude Ag = are variational functions. We use 
the parametrization 

eg = -2t[cos{qx) + cos{qy)] - 4i' cos{qx) cos{qy) - ^ 
Aq = Ai{cos{qx) - cos{qy)], 

on the 2D square lattice and the parametrization 

eg — —2ti cos(q) — 2t2 cos(2q) — /i 
Ag = Ai cos(g) A2 cos(2g) + A3 cos(3q). 



on the ID lattice. The next-nearest-neighbor hopping 
parameters, as well as the effective chemical potential fi 
and the pairing fields, are variational parameters to be 
optimized, t = ti = 1 is kept fixed to set the energy scale. 
The correlated state j^f), without backflow, is then given 

by 

1*) =^|BCS), (3) 

where 

J = exp(-i ^w.j-n.nj), (4) 

is a density-density Jastrow factor (including the on-site 
Gutzwillcr term Ui,i), with the Vij being optimized in- 
dependently for every distance \i — j\. Notably, within 
this kind of wave function, it is possible to obtain a pure 
(i.e., non-magnetic) Mott insulator by considering a suf- 
ficiently strong Jastrow factor, ^'^ i.e., Vq ~ 1/q^ {vq being 
the Fourier transform of Vij) and a metallic state when- 
ever Vq ~ 1/q. In fact, it has been shown that Vq ~ l/f?^ 
is necessary to have a vanishing quasi-particle weight in 
fermionic modelSjiS or to have a vanishing condensate 
fraction in bosonic models.— Nonetheless, much more 
difficult is to demonstrate by numerical simulations that 
such a Jastrow term is sufficient to obtain exponential 
correlations, which are suitable for a fully gapped Mott 
insulator. 

As recently shown, ^ the projected BCS state \'^) can 
be rather poor for large on-site interactions, in pres- 
ence of frustration. In order to overcome this limita- 
tion, a further improvement of the variational wave func- 
tion is needed. One possibility (which is discussed in the 
next section) is to consider the strong-coupling approach, 
based upon the S-matrix expansion .iSJ^^ Unfortunately, 
this method cannot be easily handled by quantum Monte 
Carlo simulations on large systems. However, an alter- 
native, but somewhat related, approach may be defined 
in terms of backflow correlations, which is stable and ac- 
curate even for large system sizes. 

III. THE LARGE-;/ LIMIT AND THE 
S-MATRIX EXPANSION 

Let us consider the limit of U/t — )• 00. Here, there 
are no doublons and the wave function can be written 
as l^'oc) = 'P'g|MF), where the full Gutzwillcr projector 
Vg = Y\ii^ ~ n-i^fTii^i) removes all double occupancies, 
while |MF) denotes an uncorrelated mean-fleld state, 
such as the |BCS) state introduced in the previous sec- 
tion. Then, the strong-coupling approach allows one to 
define and construct a variational wave function for the 
Hubbard model at large but finite values of U/t by 

|*5)-e-^|*oo), (5) 

where the operator S can be determined by using a 
recursive scheme. At the lowest order in t/U, iS = 
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(r+ - T-)/U, where T+ and T" correspond to the ki- 
netic terms that increase and decrease the number of dou- 
bly occupied sites by onet^^^ 

However, the wave function of Eq. ^ is hard to han- 
dle (both analytically and numerically), since S is non- 
diagonal in the basis where the electrons have defined 
positions in the lattice (in contrast to the Jastrow term, 
which is diagonal). In the largc-t/ limit, one can further 
expand the exponential and obtain: 

\^s)^{l^^S)\^oo)■ (8) 

This state has non- vanishing weights only for electronic 
configurations with zero or one doublon. Indeed, l^'cxj) 
does not contain doubly occupied sites and S may gener- 
ate at most one doublon. Therefore, after the expansion 
of the exponent, the wave function is no longer size con- 
sistent. Let us see what can be learned from this strong- 
coupling approach. Suppose that an electronic config- 
uration |a;o) has no doubly occupied sites, then, since 
'PgIxo) ~ \xo}, from Eq. (|8]) we have that: 

(xol^s) = (a:o|MF); (9) 

on the other hand, if the configuration contains one 
doublon (with the holon in one of its neighboring sites) , 
we have that: 

(x.l'fs) = ±^(2/o|MF) ± ^(zolMF), (10) 

where \yo) and |zo) are the (only) two possible configura- 
tions with no doubly occupied sites that are connected to 
|xi) by S; signs depend upon the convention used for la- 
beling electronic configurations. Therefore, the weight of 
a configuration with one doublon is related to the weight 
of two configurations without doublons (times t/U). 

This procedure for defining a variational wave function 
is exact for U ^ t, when at most one doublon is present 
in the electronic configurations; however, one is mostly 
interested into the case where an arbitrary number of 
doublons and holons are present (relevant for the generic 
case with U ^ t), and the straightforward generalization 
of this formalism becomes much more elaborated and dif- 
ficult to implement in numerical calculations. 

IV. BACKFLOW CORRELATIONS 

A. Definition of backfiow terms 

As introduced in Ref. backfiow correlations modify 
the single-particle eigenstates tpkifi.cr) of the mean-field 



Hamiltonian, like for example 'Hbcs defined in Eq. ([2]), 
according to the electronic configuration on the lattice: 

'/'^(n.a) = e0fc(n,(T) + '71 X! DiHj(j)kirj^„) 

jn.n.i 

+ in E Ai/,0fe(r„.), (11) 

jn.n.n.i 

where e = e if the site i is doubly occupied and is sur- 
rounded by at least one empty site, while e = 1 in all the 
other cases. Here e,rii,(l — 1,2) are variational param- 
eters to be optimized, Di = Ui^^rii^i and Hi = hi^-^hi,^, 
with /li.o- = 1 — rii^a. Moreover, the shorthand nota- 
tions n.n. and n.n.n. indicate nearest- and next-nearest- 
neighbor sites, respectively. In this way, already the de- 
terminant part of the wave function includes correlation 
effects, strongly improving the accuracy of the many- 
body state. When backfiow correlations are included 
in the Slater determinant, the state will be denoted by 
l^*^). This is a substantial improvement with respect to 
Jastrow factors, where electron-electron correlation is in- 
cluded only via a multiplicative term (i.e., giving a "clas- 
sical" potential acting on the electronic configuration) . 

The results for the S-matrix expansion may be com- 
pared with the ones of backfiow correlations for large 
electron-electron interactions, namely when at most one 
doublon is present. Indeed, whenever no doubly-occupied 
sites are present, we have that 

(xol*') = (xolMF) (12) 

(where we dropped the Jastrow weight, since all sites are 
singly occupied, and it represents a multiplicative con- 
stant). We discuss now the case of a single doublon and 
note that the backfiow corrections pT|) are valid for both 
up- and down-spin orbitals. We consider here, for sim- 
plicity, only nearest-neighbor backfiow terms, with pa- 
rameters e and 77 = 771 , and obtain: 

(2:1!*'') = J{xi) [±?7e(yo|MF) ±r;e(zo|MF) 

+ (xi IMF) +772(^1 IMF)] , (13) 

where, in analogy with the S'-matrix expansion, backfiow 
correlations act on configurations with doubly occupied 
sites (i.e., \xi)), transferring some weight to configura- 
tions without doubly-occupied sites (|yo) and \zq) are 
the same two configurations that appear in the strong- 
coupling expansion) ; in addition there is some weight also 
coming from configurations with one holon-doublon cou- 
ple (i.e., the original configuration |a;o) and a configura- 
tion where the two electrons of the doubly-occupied site 
have been transferred to the empty site |a;i)). The signs 
of the second and third terms of Eq. are the same as 
the ones appearing in Eq. (jlOp . The relative weights of 
the admixture in Eq. is controlled by the variational 
parameters 77 and e. 

Therefore, our way to define the backfiow correlations 
is strictly related to the strong-coupling approach; more- 
over, it allows us to consider many couples of holons 
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and doublons, with a moderate computational cost, since 
they just re-define the Slater determinant of the wave 
function. Finally, we would like to mention that the 
backflow wave function, on the contrary to the strong- 
coupling approach of Eqs. ([U, ([9]), and (fTO|) . is size con- 
sistent. 

The terms in Eq. (jlip are the dominant contributions 
to the backflow operator, having a qualitative influence 
on the properties of the correlated wave functions. In 
addition, we can also take into account further terms 
that are useful in the intermediate-coupling regime and 
correspond to all possible hopping processes: 

jn.n.i 

jn.n.n.2 

Pi ^ ni^ahi-anj-ahj^a(l)k{rj,cr)-\- 
jn.n.i 

jn.n.n.i 

jn.n.i 

jn.n.n.i 

where, in addition to e and 77/, Pi,^i{l = 1,2) are varia- 
tional parameters to be optimized. All results presented 
here are obtained by fully incorporating the backflow cor- 
rections and optimizing individually™ every variational 
parameter in and of Eq. ([2]) , in the Jastrow factor 
of Eq. (U), as well as for the backflow corrections. 
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FIG. 1: (Color online) Upper panel: Ground-state energy as 
a function of the range of backflow parameters on the 2D 
square lattice for t' /t = 0.4 and 0.75 at U/t — 8. Range 
equal to zero corresponds to no backflow correlations in the 
wave function. Lower panel: Optimized backflow parameters 
as a function of the distance |i — i|. Results for U/t — 3 
and t' /t = 0.4 correspond to a metallic state, the others to 
insulating states. Data are shown for a lattice with L = 98 
sites. 



B. Long-range correlations 

We consider now the extension of the backflow pa- 
rameters to further distances. In practice, we allow 
in Eqs. ([Tl]) and different backflow parameters 

r7|,_j| , /3|i_j| , 7|i_j| for every distance of the lattice. We 
present the results in Fig. [T] for the frustrated 2D square 
lattice. Here, it is clear that all backflow parameters are 
irrelevant for all distances larger than few lattice spac- 
ing; in particular, they are vanishing for all distances 
not connected by the hopping amplitude. Most impor- 
tantly, a remarkable gain in energy is obtained only by 
a suitable optimization of the backflow parameters up to 
second neighbors, while poor results are obtained when 
there is no backflow (or its range is smaller than the one 
of the hopping) . These results suggest that backflow cor- 
relations are indeed necessary to successfully describe the 
effect of frustrating couplings. Furthermore, we show in 
Fig. [2] the behavior of ij^i-ji for the ID lattice up to the 
fourth distance. Also in this case, the optimized back- 
flow parameters become much smaller for distances not 
connected by the hopping amplitude. 
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FIG. 2: (Color online) Optimized backflow parameters ?7|i_j| 
as a function of the distance for the ID lattice. The point 
at U/ti = 5 lies in the metallic phase, while the other points 
are located in the insulating region of the phase diagram. 
Data are shown for a lattice with L = 120 sites. 
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TABLE I: Ground-state energies for the Hubbard model on 
the square lattice with t' jt — 0.75 and L — 162, obtained 
by means of variational (VMC) and Green's Function Monte 
Carlo (GFMC) methods. Data are shown both in presence of 
and without backflow correlations. 



U/t VMC NO backflow GFMC NO backflow 



8 


-0.4048(1) 


-0.5072(3) 


10 


-0.3409(1) 


-0.4223(3) 


12 


-0.2970(1) 


-0.3617(2) 


16 


-0.2373(1) 


-0.2816(2) 


U/t 


VMC backflow 


GFMC backflow 


8 


-0.5073(1) 


-0.5307(2) 


10 


-0.4282(1) 


-0.4455(2) 


12 


-0.3704(1) 


-0.3831(2) 


16 


-0.2901(1) 


-0.2976(1) 



FIG. 3: (Color online) The energy (in unit of Atl/U), as a 
function of ti/U, for the ID Hubbard model with t2/ti = 1. 
Filled (empty) circles denote the results with (without) back- 
flow correlations. The DMRG energy of the corresponding 
Heisenberg model with J2/J1 = 1 is shown by an arrow 
Data are presented for a L = 120 lattice size. 



The effect of backflow correlations in the large U limit 
is presented in Fig. [31 where energies for large U/ti are 
shown for the frustrated ID lattice with 12/ ti = 1. Only 
the presence of backflow correlations in the wave function 
allows for a proper extrapolation to the Heisenberg limit. 
Here, the energy of the Heisenberg model is provided by 
density-matrix renormalization group (DMRG) calcula- 
tionSf2£ which are numerically exact for ID models. 

In order to assess the accuracy of backflow correla- 
tions, we compare the energies obtained within the vari- 
ational Monte Carlo (VMC) approach of Sec. |TT] with 
the data obtained by means of Green's Function Monte 
Carlo (GFMC) within the fixed-node approximation.^^ 
This approach systematically improves the variational 
wave function, extracting the best ground-state approx- 
imation with the same nodes of the starting variational 
ansatz. The data presented in Table HI for the insulat- 
ing phase of the frustrated square lattice, show that the 
variational energy, in presence of backflow correlations, 
is lower than the GFMC energy, when the starting vari- 
ational ansatz does not include the backflow term. This 
means that backflow correlations are necessary in order 
to properly describe the ground-state wave function of 
the frustrated model, capturing its correct signs. In Ta- 
ble [Hi we show also data for the unfrustrated square lat- 
tice with t' — 0. Here, backflow correlations still strongly 
improve the accuracy of the variational ansatz, but are 
less crucial in reproducing the correct signs of the ground- 
state wave functions. Indeed, for t' = 0, the signs of the 
= J7|BCS) wave function become exact in the limit 
U/t ^ Qo (e.g., the GFMC results with or without back- 
flow correlations are very close for large U/t). 



TABLE II: The same as in Table H] but with t' = 0. 



U/t VMC NO backflow GFMC NO backflow 



6 


-0.5307(1) 


-0.6444(1) 


8 


-0.3955(1) 


-0.5045(2) 


10 


-0.3353(1) 


-0.4203(3) 


12 


-0.2930(1) 


-0.3634(3) 


U/t 


VMC backflow 


GFMC backflow 


6 


-0.5961(1) 


-0.6448(3) 


8 


-0.4803(1) 


-0.5152(3) 


10 


-0.4022(1) 


-0.4278(2) 


12 


-0.3451(1) 


-0.3651(2) 



C. Effect of backflow terms on the Jastrow factor 

Here, we would like to discuss the effect of the backflow 
corrections on the behavior of the Jastrow factors across 
the metal-insulator transition. As shown in Fig. |31 the 
metallic phase is characterized by Vq ~ l/g, while the in- 
sulating phase exhibits a Vq ^ 1/?^ behavior. A further 
logarithmic divergence at g — is expected to occur in 
2Diiiii^ More details on the relation between the Jas- 
trow factor and the conduction properties may be found 
in Appendix [Xj The presence of backflow correlations 
in the wave function strongly reduces the strength of the 
Jastrow factor in the insulating phase, even if a ~ l/*?^ 
behavior is always necessary in order to induce a metal- 
insulator transition at a finite U/t. Indeed, a variational 
wave function including only a local (soft) Gutzwiller fac- 
tor P = exp(— u 7^^.t^^^4) and backfiow correlations is 
found always in the metallic phase. 



V. METAL-INSULATOR TRANSITION AND 
THE INSULATING PHASE 

Backflow correlations are a powerful tool to describe 
the whole insulating phase, from the laige-U limit down 
to the metal-insulator transition. In the following, we 
consider the density of doubly occupied sites D = 
1/L(^j ni^^rii^i), which is proportional to the interaction 
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FIG. 4: (Color online) Upper panel: Fourier transform of the 
optimized Jastrow factor Vq multiplied by \q\^ as a function 
of \q\, along the (1, 1) direction of the Brillouin zone of the 
square lattice, at t' /t = 0.75 for U/t — 7 (triangles), U/t — 
8 (squares) and U/t = 10 (circles) . Full (empty) symbols 
refer to the presence (absence) of backflow correlations, while 
different colors correspond to different lattice sizes: red (L = 
98), blue (L = 162), black (L = 242). Lower panel: The 
same quantity for a ID lattice with t2/ti = 0.9. Data are 
shown for U/ti = 5.4 (squares), U/ti = 5.6 (circles), U/ti = 
6 (up-triangles) , U /ti — 8 (down-triangles) and U /ti — 10 
(diamonds) on a L = 120 lattice size. Full (empty) symbols 
refer to the presence (absence) of backflow correlations. 



energy per site. In addition, we split the kinetic energy 
into two different parts and we compute them separately: 
/Co is the part where the hopping process does not change 
the total number of holon-doublon couples, while IC± de- 
notes the contributions where one holon-doublon couple 
is created or destroyed (see Fig. [5]). 

The results for D and the kinetic terms are presented 
in Figs. [HI [3 and [5] as a function of U for the frustrated 
square lattice with t' /t — 0.75, the unfrustrated one with 
t' — 0, and the ID lattice with t2/ti — 0.9, respectively. 
The effect of backflow can be clearly seen in the behav- 
ior of ]C±, which is systematically improved, especially 
in the insulating phase. This trend can be understood 
by looking at Eq. (ITT|) : as long as a holon-doublon cou- 
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FIG. 5: (Color online) Upper raw: Hopping terms of a spin 
up electron from site i to site j that contribute to ICo- Lower 
raw: hopping terms of a spin up electron from site i to site j 
that contribute to /C±. Lattice sites i and j are represented 
as boxes and spin up/down electrons inside them represent 
the electronic configurations. 



pie is formed, backflow correlations favor its recombina- 
tion into single-occupied sites, increasing the part of the 
kinetic energy that changes the total number of holon- 
doublon couples. Moreover, the effect of backflow corre- 
lations is particularly strong close to the metal-insulator 
transition, where the increasing of the holon-doublon re- 
combination washes out the jump in JC±. On the con- 
trary, the behavior of JCq shows a suppression in presence 
of backflow correlations, since hopping terms that does 
not change the total number of holon-doublon couples 
contribute less to the total energy. We want to point out 
that, in presence of frustration, JCq exhibits a clear jump 
at the metal-insulator transition, similarly to the dou- 
ble occupation D, both for the square (see Fig. |6]) and 
the ID lattice (see Fig. HI). Instead, for the unfrustrated 
square lattice, the jump vanishes and the metal-insulator 
transition is characterized only by an inflection point (see 
Fig.IZl). 

Finally, we would like to discuss how the variational 
method allows us to estimate the charge gap in the in- 
sulating phase, by assessing the particle-hole excitations 
over the variational state. This is possible by consid- 
ering just ground-state expectation values, without di- 
rectly calculating energy differences. Indeed, within the 
context of the single-mode approximation (SMA), which 
was proposed for the liquid Helium^^ and further ap- 
plied to fermionic systems^^i^i it is possible to find 
a relation between the particle-hole excitation energy 
and the static structure factor N(q) — (n^qUq) where 
Uq = ^e^'^^nr.cr IS the Fourier transform of the 

particle density. Indeed, a variational ansatz of the low- 
est energy state with a given momentum q, can be 
obtained by applying Uq to the ground state wave func- 
tion 1 5") (or one approximation for it), namely: 



\^,)=nq\^) 



The charge gap for the limit q 
Appendix iBl and results: 



(15) 

is then derived in 



Eq = -^(lur.ML 
' 2d \q^a N{q) 



[(Ari)2/Ci + (Ar2)2/C2] , (16) 
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FIG. 6: (Color online) Upper panel: Density of doubly occu- 
pied sites D as a function of U jt. Lower panel: /Co (circles) 
and K,± (squares) as a function of U jt. Full (empty) sym- 
bols refer to the presence (absence) of backflow correlations in 
the wave function. The metal-insulator transition is marked 
by an arrow. Data refer to a frustrated square lattice with 
i! It = 0.75 and L = 162. 

where d is the dimensionality, K\ and /C2 are the nearest- 
and the next-nearest-neighbor kinetic energy per site, 
while Ari and Ar2 are the distances between nearest- 
and next- nearest- neighbors on the lattice (e.g., Ari = 1 
and Ar2 = 2 or \/2 on the ID or 2D square lattice). 

The metallic phase is characterized by Af((z) ~ q for 
q — >■ 0, which implies a vanishing gap for particle-hole 
excitations. On the contrary, in the insulating phase, 
iV(q) ^ for g — > 0, enclosing the fact that the charge 
gap is finite. 

In Fig. ini we show the gap for g ^ for the ID Hub- 
bard model with ^2/^1 = 0.9. The gap closes at the 
metal-insulator transition, as expected, and grows up lin- 
early with U for large values of U jti. For comparison, 
we also show the same gap in the case <2 = (where a 
Mott insulating state takes place for infinitesimal values 
of C//ii). 



VI. CONCLUSIONS 



FIG. 7: (Color online) The same as in Fig. [1 but for t' = 0. 
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The search for high-quality variational descriptions of 
the ground-state properties of correlated and frustrated 
electron systems is a long-standing issue. The problem 



FIG. 8: (Color onhne) The same as in Fig. [6l but for a ID 
lattice with t2lt\ = 0.9 and L = 120. 
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FIG. 9: (Color online) Charge gap Eg for g -> for the ID 
Hubbard model with i2/ti = 0.9, on L = 120 and L = 240 
lattices. Results without backflow correlations for the L = 
120 case are shown for comparison. Inset: charge gap for the 
ID Hubbard model with t2 = on a L = 120 lattice size. 
Axis labels are the same of the main plot. 



lies in the correct description of emergent energy scales 
of order t'^/U, which are responsible both for the gen- 
eration of antiferromagnetic correlations and spin-liquid 
behaviors of frustrated systems, as well as for incipient 
and real superconducting state formation. 

The description of "dynamical" energy scales, viz of en- 
ergy scales which are not explicitly present in the original 
Hamiltonian but are generated by dynamical processes 
(e.g., the super-exchange J = At^ /U), is a challenge for 
both variational approaches and, more generally, for nu- 
merical simulations. It can be circumvented in the limit 
of large interactions, i.e., ?7 3> i, by using the S'-matrix 
expansion outlined in Sec. IIIII One can then transform 
the Hubbard model into the t—J model, for which there 
is no need to generate the antiferromagnetic energy scale, 
since it is already present explicitly in the model. 

The t—J model approach has played a key role in devel- 
oping the RVB approach to high-temperature supercon- 
ductivity. Within this framework, it is however not pos- 
sible to study and correctly describe the Mott-Hubbard 
transition. With the backflow corrections, investigated in 
detail in the present study, it is instead now possible to 
describe, with a high degree of accuracy, the dynamical 
energy-scales of the Hubbard model both in the large- 
and the intermediate- J7 region. We hence believe that 
substantial progress has been made towards the solution 
of a long-standing issue, namely the correct description 
of emergent energy scales in correlated and frustrated 
electron systems. 

L.F.T. and C.G. acknowledge the support of the Ger- 
man Science Foundation through the Transregio 49. 



Appendix A: Jastrow factors and charge gap 

As discussed in Sec. [Vj the insulating or metallic na- 
ture of the ground state of the Hubbard model can 
be determined by looking at the static structure factor 
N{q) = {ri-qUq). Indeed, the behavior of this quantity 
for (7 — is related to the presence of a gap in the charge 
excitations of the system by Eq. (IT51) . The insulating 
phase is characterized by N{q) ~ for g — 0, while in 
the metallic phase N{q) ~ q. 

A relation between the static structure factor and the 
Jastrow factor = exp(— 1/2 VqUqU^q) has been de- 
rived in the context of liquid Helium by Reatto and 
Chester,^? using a Gaussian approximation for the proba- 
bility density associated to the uncorrelatcd ground state 
wave function |MF). We have that: 



N{q) = 



N\q) 



l + 2vqm{qy 



(Al) 



where N^{q) is the static structure factor for the uncor- 
related wave function |MF). This relation is rigorously 
valid in the weak-coupling regime; nevertheless, a similar 
behavior is found in the (strong-coupling) ID insulat- 
ing phase. The optimization of Vi,j for every distance 
|z — j|, will lead, in g-space, to Vq ^ \/q for the metalhc 
phase (in any dimension) and to Vq ~ l/?^ (in ID) and 
Vq ^ \og{q)/q^ (in 2D), see Fig. g] Since N^{q) ^ const 
for |MF) = |BCS) and Vq diverges for g — > 0, the con- 
dition 2vqN^{q) ^ 1 holds in the limit g — > and one 
obtains 



N{q) ^ - 
v„ 



(A2) 



Therefore, the small-g behavior of the structure factor 
reflects the smaW-q behavior of the Jastrow factor. In the 
2D insulating phase, although Vq ^ \og{q)/q^ is found, 
the static structure factor shows a quadratic behavior, 
i.e., N{q) ~ q^ , highlighting the fact that corrections to 
Eq. (|A1[) are present in the 2D Mott phase. 



Appendix B: Derivation of the single-particle charge 

gap 

Given the ansatz for the excited state of Eq. (|15p . 



where Un 



is the Fourier- 



transformed particle density, the variational estimator of 
the excitation energy is then given by 



{^q\{H-E^)\^q 



mn^q[H,nqm 



(Bl) 



where T-L is the Hamiltonian of Eq. ([T]) . The sum of both 
commutators (n-qHnq — n^qUqH) + {n^qHnq ~ Hn^qUq) 
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is equivalent to the double commutator 



Therefore, in the limit of q 0, Eg can be written as: 



[n^q, [H,n,]] = n^q{T-Lnq - UqH) - {Unq - nq'H)n^q 

(B2) 

due to the inversion symmetry q -H- —q. Therefore, the 
excitation energy is given by 



E„ 



1 {nn-„m,nq]m 

2 Nn 



(B3) 



where N{q) = {"i'\n_qnq\'^) is the static structure fac- 
tor for the ground state. The double commutator 
[H, Hq]] can be straightforwardly evaluated and in- 
volves the kinetic term only, since the potential term 
of the Hamiltonian contains density-density interactions 
that commute with n^: 

[n^q, [H, nq]] = J ^(efe+, + e^-q - 2efe)4^^c,^^. (B4) 



E„ 



1 fy 

lim 

2d ^9^0 N{q) 



((Ari)2/Ci + (Ar2)2/C2) 



(B5) 

where d is the dimensionality, ICi and /C2 are the nearest- 
and the next-nearest-neighbor kinetic energy per site, 
while Ari and Ar2 are the distances between nearest- 
and next-nearest-neighbors on the lattice (e.g., Ari = 1 
and Ar2 = 2 or a/2 on the ID or 2D square lattice). 
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